home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
The World of Computer Software
/
The World of Computer Software.iso
/
yamp16.zip
/
TESTGRAF.CPP
< prev
next >
Wrap
C/C++ Source or Header
|
1993-01-11
|
6KB
|
208 lines
/*
YAMP - Yet Another Matrix Program
Version: 1.6
Author: Mark Von Tress, Ph.D.
Date: 01/11/93
Copyright(c) Mark Von Tress 1993
Portions of this code are (c) 1991 by Allen I. Holub and are used by
permission
DISCLAIMER: THIS PROGRAM IS PROVIDED AS IS, WITHOUT ANY
WARRANTY, EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED
TO FITNESS FOR A PARTICULAR PURPOSE. THE AUTHOR DISCLAIMS
ALL LIABILITY FOR DIRECT OR CONSEQUENTIAL DAMAGES RESULTING
FROM USE OF THIS PROGRAM.
*/
#include "virt.h"
//required global declaration for the
// matrix stack object
// unsigned int _stklen = STACKLENGTH;
MStack *Dispatch = new MStack;
VMatrix &getx(int N)
// create an x matrix
{
Dispatch->Inclevel();
VMatrix x, c1, x2;
c1 = Fill(N, 1, 1.0);
x = Index(1, N) - ((double) N) *0.5;
x = Ch(c1, x);
// push x onto the stack
Dispatch->Push(x);
// decrement the subroutine nesting level
// and return the stack top
return Dispatch->DecReturn();
}
VMatrix &gety(VMatrix &x, VMatrix &beta)
// create a y matrix
{
Dispatch->Inclevel();
VMatrix y;
y = x*beta;
srand(123);
for (int i = 1; i <= y.r; i++) {
// use sum of 3 uniforms for an approximate
// normal random variable
int u = random(100) + random(100) + random(100) + 3;
y.M(i, 1) = y.m(i, 1) + 5.0*sin(3.14*((double) (i % 8)) / 7.0)
+ ((double) (u - 150)) / 300.0;
}
Dispatch->Push(y);
return Dispatch->DecReturn();
}
VMatrix ®ression(VMatrix& x, VMatrix& y)
// do a multiple linear regression
{
Dispatch->Inclevel();
VMatrix yx, reg, betahat;
int N = x.r, p = x.c;
// solve for regression parameters using sweep
yx = Ch(y, x);
reg = Sweep(2, p + 1, Tran(yx) *yx);
// calculate mean square error
reg.M(1, 1) = reg.m(1, 1) / ((double) (N - p));
reg.DisplayMat();
// solve regression using normal equations
betahat = Inv(Tran(x) *x) *Tran(x) *y;
Dispatch->Push(betahat);
return Dispatch->DecReturn();
}
void plotResiduals(VMatrix &resids)
{
Dispatch->Inclevel();
VMatrix grf = Ch(Index(1, resids.r), resids);
GMatrix Agraph(grf, - '%');
*Agraph.PathToDriver = "C:\\tc\\bgi";
*Agraph.title = "Residuals for data";
*Agraph.title2 = "with serial correlations with frequency 0.25";
*Agraph.yname = "Residuals";
*Agraph.xname = "Index";
Agraph.Href(0.0);
Agraph.Show();
Dispatch->Cleanstack();
Dispatch->Declevel();
}
VMatrix &GetSerialCovar(VMatrix &R, VMatrix &spectdens)
{
// Parameters to CORR in Timeslab
// correlogram = CORR(x=R,n=R.r,M=R.r-1,Q=2*R.r,
// iopt=1,r0=r0, per = spectdens)
// covar = correlogram*r0
Dispatch->Inclevel();
VMatrix centered, z, covar;
double n = (double) R.r;
// center a column vector
centered = R - Sum(R).m(1, 1) / n;
// zero pad to length 2n: 2n periodic for full
// sample spectral density
centered = Cv(centered, Fill(R.r, R.c, 0));
// take fft
z = Fft(centered);
// take convolution : gives sample spectral density
spectdens = Sum(z % z / n, COLUMNS);
// inverse fft for serial correlation (autocorrelation)
covar = Fft(spectdens, FFALSE);
// throw away last half.
covar = Submat(covar, R.r, 2);
Dispatch->Push(covar);
return Dispatch->DecReturn();
}
void plotCorrelogram(VMatrix &serial)
{
Dispatch->Inclevel();
int n = serial.r;
double sigma = serial.m(1, 1);
VMatrix Correlogram = serial / sigma;
Correlogram = Submat(Correlogram, n, 1, 2, 1);
VMatrix graf = Ch(Index(1, Correlogram.r), Correlogram);
GMatrix Agraph(graf);
*Agraph.PathToDriver = "C:\\tc\\bgi";
*Agraph.title = "Serial Correlations";
*Agraph.title2 = "for sample residuals";
*Agraph.yname = "Serial correlations";
*Agraph.xname = "Lags";
Agraph.Href(0.0);
Agraph.Show();
Dispatch->Cleanstack();
Dispatch->Declevel();
}
void plotPeriodogram(VMatrix &spectdens)
{
// calculate a standardized periodogram on the log scale
Dispatch->Inclevel();
int n = spectdens.r;
// this works because data is already centered, which
// forces spectdens.m(1,1) = 0.0;
double sigma = Sum(spectdens).m(1, 1) / ((double) n);
VMatrix Periodogram = spectdens / sigma;
Periodogram = Mlog(Submat(Periodogram, n / 2 + 1, 1, 2, 1));
// frequencies
double dn =(double) (n / 2 + 1);
VMatrix graf = Ch(Index(1, Periodogram.r) / dn, Periodogram);
GMatrix Agraph(graf);
*Agraph.PathToDriver = "C:\\tc\\bgi";
*Agraph.title = "Periodogram";
*Agraph.title2 = "for sample residuals";
*Agraph.yname = "Log periodogram";
*Agraph.xname = "Frequencies";
Agraph.Vref(0.25);
Agraph.Show();
Dispatch->Cleanstack();
Dispatch->Declevel();
}
main()
{
Dispatch->Inclevel();
VMatrix x, beta("beta", 2, 1), y, betahat, resids, serial;
Setwid(15);
Setdec(10);
beta.M(1, 1) = 1;
beta.M(2, 1) = 0.5;
x = getx(128);
y = gety(x, beta);
betahat = regression(x, y);
betahat.Nameit("Text book betahat");
(Tran(beta)).DisplayMat();
(Tran(betahat)).DisplayMat();
resids = y - x*betahat;
plotResiduals(resids);
VMatrix spectdens;
serial = GetSerialCovar(resids, spectdens);
plotCorrelogram(serial);
plotPeriodogram(spectdens);
vclose();
}